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Relaxation dynamics of complex quantum systems with strong interactions towards the steady state is a 
fundamental problem in statistical mechanics. The steady state of subsystems weakly interacting with their 
environment is described by the canonical ensemble which assumes the probability distribution for energy to 
be of the Boltzmann form. The emergence of this probability distribution is ensured by the detailed balance 
of the transitions induced by the interaction with the environment. Here we consider relaxation of an open 
correlated quantum system brought into contact with a reservoir in the vacuum state. We refer to such a system 
as emissive since particles irreversibly evaporate into the vacuum. The steady state of the system is a statistical 
mixture of the stable eigenstates arising due to the binding energy. We found that, despite the absence of the 
detailed balance, the stationary probability distribution over these eigenstates is of the Boltzmann form in each 
A-particle sector. A quantum statistical ensemble corresponding to the steady state is characterized by different 
temperatures in the different sectors, in a contrast to the Gibbs ensemble. We investigate the transition rates 
between the eigenstates to understand the emergence of the Boltzmann distribution and find their exponential 
dependence on the transition energy. We argue that this property of transition rates is generic for a wide class of 
emissive quantum many-body systems. 

PACS numbers: 05.30.Ch, 05.70.Ln 


Whether and how a quantum system brought out of equilib¬ 
rium reaches its steady state is a fundamental question which 
has recently attracted much attention KM. Expectation val¬ 
ues of local observables after relaxation are typically deter¬ 
mined by integrals of motion while the memory of micro¬ 
scopic details of an initial state is lost. In quantum systems the 
mechanism of thermalization is rooted in the properties of in¬ 
dividual eigenstates as stated in the eigenstate thermalization 
hypothesis (ETH) ITHD. The basic statement of ETH for non- 
integrable systems is the smoothness of eigenstate expectation 
values of local observables as functions of eigenenergies. Eor 
systems with integrals of motion this statement remains true if 
eigenstates are taken from the same subspace IfTOin^ . In the 
steady state the density matrix of any subsystem is diagonal 
and its elements are the same for the typical initial eigenstates 
of the full system iBHia- ETH may break down for rare 
eigenstates in the low-energy part of the spectrum iniiiii. 

A subsystem weakly coupled to the rest of the system can 
be viewed as an open system with its surroundings acting as a 
reservoir. Transitions induced by the coupling to the reservoir 
govern relaxation dynamics of the open system to the steady 
state which is described by the (generalized) Gibbs ensemble 
EMa. These transitions satisfy the detailed balance prin¬ 
ciple in the steady state and the probability distribution over 
many-body eigenstates of the open system is of the Boltzmann 
form. However, this general scenario does not describe a spe¬ 
cial case of an open correlated quantum system brought into 
contact with a reservoir in the vacuum state. In this case par¬ 
ticles evaporating from the system never return, and therefore 
the detailed balance principle does not hold. In the follow¬ 
ing we refer to such a system as emissive. Although particles 
can only leave the system, its steady state may be populated if 
there is a binding energy for escaping particles. Such a pop¬ 
ulated steady state is a statistical mixture of multiple eigen¬ 
states whose probability distribution is not a priori known. 



Figure 1: The stationary probability distribution of the emissive 
system of hard-core bosons on the 14-site lattice (geometry is shown 
in the inset) is shown. The initial state of the system is the max¬ 
imally occupied pure state with No = 14. In sectors with at least 
3 stable states the data is fitted with the Boltzmann distribution 
Pf oc exp(-E„/Tfj). The temperatures Tjv are presented in the ta¬ 
ble. 


The quantum statistical ensemble for emissive systems de¬ 
fined in this way has never been considered so far. 

Open systems coupled to a vacuum reservoir appear ubiq¬ 
uitously in various fields of natural sciences including surface 
science ll^l24ll . quantum optics Il25ll26l . nuclear physics EtII 
and astrophysics ESll . In the field of cold atoms an established 
experimental example is a trapped atomic or molecular gas in 
a vacuum chamber mM- Collisions of trapped particles re¬ 
sult in internal equilibrium of the gas which alters as particles 
escape. The loss of particles can be accompanied by cool¬ 
ing. This process of evaporative cooling is a key technologi¬ 
cal development used to achieve Bose-Einstein condensation 
of cold atoms in magneto-optical traps Il32l435l . The systems 
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studied in experiments with ultracold atoms typically contain 
from 10^ to 10® particles. Statistical properties of these sys¬ 
tems are usually probed by measuring their local observables. 
Recent experimental progress has made it possible to study 
smaller systems Il36l - l39l . For these systems a direct measure¬ 
ment of the probability distribution over many-body states can 
be realized. In particular, it would allow to probe statistical 
properties of the steady state of emissive quantum systems. 

In this paper we study an emissive quantum system of hard¬ 
core bosons on a lattice. We demonstrate that the probability 
distribution in the steady state is of the Boltzmann form in 
each -particle sector. We characterize this steady state by 
a quantum statistical ensemble and discuss its application to 
calculating expectation values of observables. We show that 
the physical mechanism behind the emergence of the Boltz¬ 
mann distribution is rooted in the behaviour of the transition 
rate. This rate appears to be a regular function of the transi¬ 
tion energy though there is no detailed balance with the reser¬ 
voir. We connect this behaviour with statistical properties of 
off-diagonal matrix elements of local annihilation operators 
entering expressions for transition rates. 

We study a system of hard-core bosons on a two- 
dimensional lattice coupled to the vacuum reservoir with the 
Hamiltonian H - Hg + Hr + Hj. The lattice is described by 
the Hamiltonian 

Hs^-Y,h,j{b\bj + b]bi) (1) 

(ij) 

with a constraint that each site can be occupied with no more 
than one particle. Here b^. (bi) is a creation (annihilation) op¬ 
erator on the site i, hij are hopping amplitudes and the sum 
runs over the pairs of the nearest-neighbour sites. For calcu¬ 
lations we use a lattice of 14 sites (see inset of Figure 1). To 
ensure non-integrablility of the lattice, the different values are 
assigned to the hopping amplitudes hij. These values are taken 
from the interval between 0.8 and 1.2 (in arbitrary units of en¬ 
ergy). The vacuum reservoir is described by the Hamiltonian 
Hr — Yik with operators (a^) creating (annihilating) 

particles in reservoir modes k with energy ek- The interaction 
between the lattice and the reservoir is introduced through the 
Hamiltonian H/ - aYjkAo.\bi + b\ak)0{ek - eq), where the 
unit step function 6(sk - eq) accounts for the binding energy 
sq. This binding energy is a tunable parameter in our calcu¬ 
lations. We consider the weak coupling limit a « 1 and as¬ 
sume that escaping particles immediately lose the coherence 
with the lattice and do not return back to it. 

We describe the reduced dynamics of the lattice by a master 
equation 14011411 for its density matrix ps . The density matrix 
is represented in the basis of eigenstates of Hg . These eigen¬ 
states |A^, n) are divided into sectors according to the number 
of particles N in the lattice. Corresponding eigenenergies are 
denoted by E^. We perform leading order calculations of the 
transitions between the eigenstates generated by the coupling 
to the reservoir l42ll . In the absence of special symmetries 
of the lattice there are two constraints on these transitions: 
(i) AN = -1 which is due to the type of coupling, and (ii) 


AE < -So which is due to the binding energy for escaping 
particles. 

We consider the dynamics of the system at the coarse¬ 
grained time scale dt ~ tio^s, where is the characteris¬ 
tic time of the particle loss process. In the weak coupling 
regime tioss is much larger than the dephasing time, 

so that only the diagonal elements of the density matrix con¬ 
tribute to dynamics. These elements (A^, n|ps|A^, n) = may 
be viewed as the probability distribution of the lattice over its 
eigenstates. The master equation is reduced to the system of 
the rate equations for this probability distribution 

pV _ \ ’ pA'+lpV+l _ \ ’ nN pN /p\ 

~ / i ^nm / ; ‘^mn^n ’ 

m m 

where is the transition rate from the state |A^, n) into the 
state |A^ - l,m). For transitions allowed by the selection rules 
the rates are estimated using Fermi’s golden rule as 

L 

= 2nQ^a^ niMN, n)\^. (3) 

i=\ 

Here the step-like density of states of the reservoir Q. - 
Qo6(e - Eo) has been used that corresponds to the two- 
dimensional motion of escaping particles. The particle loss 
process can be viewed as a sequence of transitions which 
brings the system from the initial to one of the stable states 
from which no further transitions can occur (see Supplemen¬ 
tary). There is the trivial stable state with no particles and 
occupied stable states which appear due to the constraint on 
the transition energy. For calculations we use the maximally 
occupied pure state with No - 14 as the initial one and choose 
So = 0 to maximize the number of the stable states which can 
be achieved (see Supplementary). For the 14-site lattice this 
number is 26 and the maximal occupation of the stable states 
is 7. 

The probability distribution of the system in a steady state 
is shown in Figure 1. In each A^-particle sector of the Hilbert 
space (having at least 3 stable states) this stationary distribu¬ 
tion is found to be of the Boltzmann form. Parameters of the 
exponential fit of the data Ea 

Zrt 

define a quantum statistical ensemble. This ensemble is char¬ 
acterized by the set of temperatures T^ and the probabilities 
Pff for the system to have N particles in the steady state. The 
partial partition functions are summed over 

achievable stable states. We assign zero probabilities to the 
unstable states which are not present in the ensemble. Sta¬ 
tionary expectation values of the observables are given by a 
standard expression (A) = YiNnPn{N,n\A\N,n). We note that 
the steady state can be characterized by this statistical ensem¬ 
ble for an arbitrary initial state of the system (see Supplemen¬ 
tary). Variation of the initial state only changes parameters Pa? 
and Tat of the distribution. 
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Figure 2: The probability distribution of the intermediate states in 
each A^-particle sector is shown. The initial state of the system is the 
maximally occupied pure state with No = 14. For illustrative pur¬ 
poses the rates are multiplied by the factor 10'''. The stable states are 
indicated with empty circles. In each sector the points accumulating 
80% of the total probability are fitted with the exponential functions 
Pn{E) oc exp(-£/rA(). The temperatures are presented in the 
table. 

Emergence of the distribution. The emergence of the 
Boltzmann distribution in the emissive system can be under¬ 
stood from statistics of the transitions which the system un¬ 
dergoes as particles evaporate. For each sector we calculate 
the occurrence probabilities of the intermediate states, i.e. the 
states through which the system evolves from the initial to one 
of the stable states. These distributions in sectors (referred to 
as intermediate) are shown in Figure 2 as functions of the en¬ 
ergy of the states. In sectors with = 3 11 the dependence 

is smooth and fluctuations of the data are less than the size 
of points. In each sector we consider the part of the spectrum 
accumulating the probability of O.SPat and approximate corre¬ 
sponding points by the Boltzmann distribution 0. For < 7 
the data includes stable states and the temperatures are close 
to those in Figure 1. We note that measuring the occurence 
probabilities of the intermediate states will require collecting 
statistics from an ensemble of emissive systems. Assuming 
that the spectrum of the system is discrete and not degenerate, 
a sequence of the intermediate states can be directly deter¬ 
mined by measuring the energies of escaping particles. 

The smooth energy dependence of intermediate probability 
distributions can be understood from the statistical properties 
of the transition rates in the system. Figure 3 shows plot¬ 
ted as a function of the transition energy e - En - E^. The 
data for R^{s) in each sector can be approximated by a reg¬ 
ular dependence with state-to-state fluctuations. For an arbi¬ 
trary initial A-particle ensemble the probability distribution in 
(N- l)-particle sector after the loss of a particle is determined 
by 

n 


Possible state-to-state fluctuations of are convolved with 
the regular function of the transition energy and almost do not 
translate into the probability distribution in the following sec¬ 
tor. The probability distribution becomes a smooth function 
of energy after the loss of several particles if transition rates 
in neighboring sectors are uncorrelated. 

The analysis of lattices of different sizes (see Supplemen¬ 
tary) shows that the transition rate is always a regular function 
of the transition energy with state-to-state fluctuations. The 
fluctuations decrease with the increase of the size of the lat¬ 
tice. We expect that they vanish in the thermodynamic limit, 
though it is not evident from the available range of lattice 
sizes. This assumption is supported by the analysis of the 
transition rates in the thermodynamic limit presented below. 

Thermodynamic limit. Let us show that the exponential 
dependence of the transition rate on the transition energy is a 
universal feature of quantum many-body systems coupled to 
the vacuum reservoir. For this purpose we consider evapora¬ 
tion in a generic system of A » 1 particles with well defined 
energy E. Temperature T and chemical potential yu of the sys¬ 
tem are related to the density of many-body states A(N,E) 
via well known thermodynamic relations 1/T - ^InA, 
jj. - -T ^InA. It is reasonable to assume that the energy 
e of a particle escaping to the vacuum obeys the Boltzmann 
distribution with some temperature T' (not necessarily equal 
to T, e.g. because of the adiabatic cooling during evapora¬ 
tion). The probability density of the energy of the leaving 
particle is p{s) oc a(e)e“®/^', where a{e) is the density of 
single particle states in the vacuum reservoir. On the other 
hand, the same probability can be estimated with the Fermi’s 
golden rule as p{s) oc R^(s)A(N - 1,E - s). One can use 
the first order Taylor series expansion of In A and the ther¬ 
modynamic relations above to show that A{N - 1,E - e) = 
A(N,E)e~^‘^~^^^^. Comparing two expressions for p{s), we 
conclude that Inf?(e) = -{1/T' - 1/T)s-\na{s). For our spe¬ 
cific case of two-dimensional reservoir with a{s > 0) = const 
the transition rate is a regular exponential function of the tran¬ 
sition energy, in accordance with results shown in Figure 
The difference between the inverse temperatures of the gas 
and the leaving particle is determined by the slope angle of 
In R{s) function. We conclude that an exponential dependence 
of the transition rate on the transition energy is a property of 
the wide class of emissive quantum many-body systems where 
escaping particles satisfy the Boltzmann distribution. 

The regular dependence of the transition rate on the transi¬ 
tion energy relates to statistical properties of the off-diagonal 
matrix elements {N - l,m\bi\N, n) which enter the expression 
It is important to note that A,- are non-Hermitian anni¬ 
hilation operators and the matrix elements are calculated be¬ 
tween states in the different sectors, in contrast to known stud¬ 
ies ms El. We argue that the smooth dependence of these 
matrix elements on the transition energy is similar to that of 
eigenstate expectation values in ETH. Fluctuations of the ex¬ 
pectation values vanish as the size of the system increases 
iMii. From our scaling analysis we expect that fluctua¬ 
tions of {N - \,m\bi\N, n) have the same property. 
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Figure 3: The transition rates from the states in the A^-particle 
sector to the states in the (N - l)-particle sector as a function of the 
transition energy s = - £■"“* for iV = 3,..., 12. The allowed 

transitions correspond to s > Sq = 0. For illustrative purposes 900 
points are shown for each N and the rates are multipled by the factor 
10 ". 


We note that an experimental verification of our findings 
will require measuring the probability distribution of a sys¬ 
tem under study over its many-body states ll48l . This type of 
measurements would extend the method of determining tem¬ 
perature from average values of (local) observables. While 
the latter is available both for finite-size and macroscopic sys¬ 
tems, the former is only realizable for the systems of finite 
size. 

Cooling/heating effect. Emissive systems are typically 
cooled down by the particle loss process since leaving par¬ 
ticles accumulate on average more energy than the remaining 
ones. Preparation of the system in a Gibbs state with No < 14 
allows us to study the cooling/heating effect in the system and 
its dependence on the initial temperature Tq. Fig. shows 
the inverse temperature l/T as a function of N for the Gibbs 
states with = 12 and different initial temperatures. The en¬ 
ergy spectrum of the system is bounded both from below and 
above, so the initial states with negative Tq can also be consid¬ 
ered. In the case of a negative initial temperature the system 
always cools down as particles escape. In the case of a small 
positive initial temperature the loss of the first particle leads to 
the abrupt heating of the system. This can be explained by the 
shift between the spectra of the neighbouring sectors which 
leads to the possibility of transitions from the lowest energy 
state in the Wparticle sector to the excited {N - l)-particle 
states. For the system under study this condition can be sat¬ 
isfied only for N > 7 so that no heating is observed below 
half-filling. 

Conclusions. To summarize, we demonstrate an emer¬ 
gence of the Boltzmann distribution in each A^-particle sector 
of an emissive system of hard-core bosons on a lattice. The 
steady state of the system is described by a quantum statisti- 
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Figure 4: The plot shows the dependence of the inverse temperatures 
l/Tjv on the number of particles for different initial conditions which 
are defined by the set of initial inverse temperatures 1/Tq (shown in 
the inset table). The plot demonstrates that system may undergo both 
heating and cooling process during the evaporation depending on the 
initial temperature of the system and the filling of the system. 


cal ensemble characterized by a set of temperatures. Smooth 
probability distributions are the consequence of the statistical 
properties of the transition rates in the system. The transi¬ 
tion rates are expressed in terms of the off-diagonal matrix 
elements of local annihilation operators and depend smoothly 
on the transition energy. This allowed us to draw the analogy 
with ETH which asserts the smooth dependence of the eigen¬ 
state expectation values on the energy. We expect that this fea¬ 
ture of transition rates is a generic property of the correlated 
emissive quantum systems and the formation of the smooth 
Boltzmann distribution can be observed in contemporary ex¬ 
periments with ultracold atoms. The physics demonstrated in 
this paper is strongly determined by a particular choice of the 
reservoir which can only absorb particles. This excludes ther- 
malization within the same sector due to virtual processes and 
makes low-energy populated states stable. 
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Hilbert space and transitions 


The Hilbert space of the system can be represented as the 


graph shown in Fig. S1 The vertices of this graph represent 
the many-body eigenstates \N, n) which are classified by the 
number of particles N and the energy Directed edges con¬ 
nect these vertices according to the transitions between the 
eigenstates generated by the coupling to the vacuum reser¬ 
voir. The dynamics of the system can be considered as a 
walk through the graph. Transitions should satisfy two con¬ 
straints; (i) AA = -1, (ii) AE < —eq. By possible transitions 
we mean all the transitions satisfying the first constraint. We 
divide them into allowed transitions which also satisfy the sec¬ 
ond constraint and forbidden transitions which do not satisfy 
the second constraint. We also introduce two categories for 
the eigenstates of the system. For a given initial state we de¬ 
fine an achievable eigenstate as an eigenstate which can be 
reached from the initial one through the allowed transitions. 
We define stable eigenstates as the eigenstates from which no 
further transitions are allowed. Only achievable stable eigen¬ 
states are present in the steady state of the system. We note 
that the dynamics of the system clearly does not satisfy the 
detailed balance principle since all possible transitions are ir¬ 
reversible. 



Figure SI: Graphical representation of the Hilbert space. The 

system of hard-core bosons on a lattice of 4 sites is considered. Bars 
representing many-body eigenstates \N, n) are placed into the coordi¬ 
nate system ’number of particles - energy’. Initial maximally occu¬ 
pied pure state is indicated by the green colour. Allowed transitions 
and some forbidden transitions for eo = 0 are shown with solid and 
dashed arrows correspondingly. The achievable eigenstates are in¬ 
dicated with the blue colour. The stable eigenstates are represented 
with empty bars. After the loss of particles the systems in the en¬ 
semble can be in either the lowest energy 2-particle eigenstate or the 
lowest energy 3-particle eigenstate. 


Binding energy 


The binding energy sq determines the allowed transitions 
and the eigenstates which are present in the steady state of the 
system. The dependence of the number of these eigenstates 
on Eo for the initial maximally occupied pure state is shown 


in Fig. S2 For the extreme values of the binding energy there 
is only one achievable stable eigenstate: (i) the empty eigen¬ 
state for fio = (all possible transitions are allowed), (ii) the 
maximally occupied eigenstate for eq - oo (all possible tran¬ 
sitions are forbidden). For calculations we set sq = 0 which 
corresponds to the maximal value of the number of eigenstates 
present in the steady state of the system. 



Binding energy 


Figure S2: Number of achievable stable eigenstates. The system 
of hard-core bosons on a 14-site lattice is prepared in the maximally 
occupied pure state. The plot shows the number of achievable stable 
eigenstates as the function of the binding energy. 
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Energy 


Figure S3: Stationary distributions. The system is prepared in pure 
states with No = 12 particles and different energies: (a) Eq = -3.00, 
(b) Eo = -0.06, (c) Eo = 0.06, (d) Eq = 3.00. The stationary dis¬ 
tributions are of the Boltzmann form in each Ai-particle sector for all 
these initial conditions (temperatures are presented in tables). 


Dependence on the initial energy 


As stated in the paper, the steady state of the system is char¬ 
acterized by the Boltzmann distribution in each sector. Here 
we demonstrate that the variation of the initial state of the sys¬ 
tem only changes parameters of the stationary distribution but 
not its shape. Fig. S3 shows stationary distributions for the 
initial pure states with Nq = 12 particles and the different en¬ 
ergies Eq. We observe that (i) temperatures increase as Eq 
increases, (ii) for the initial states with close energies station¬ 
ary distributions are similar. 


Scaling analysis 


The amplitude of the state-to-state fluctuations of the tran¬ 
sition rate as the function of the transition energy depends on 
the size of the system. Fig. S4 shows that the standard devi¬ 
ation of the state-to-state fluctuations becomes smaller as the 
size of the system increases. Lattices of the different sizes are 
obtained by taking out sites from the 14-site lattice. 



Figure S4: Scaling analysis. The plot shows the standard deviation 
of logio/? in the interval of transition energies [-0.1,0.1] as the func¬ 
tion of the number of sites in the lattice. Geometries of the lattices of 
different sizes are shown in the inset. 

















































